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A new third-order Energy Stable Weighted Essentially Aon Oscillatory (ESWENO) finite 
difference scheme for scalar and vector linear hyperbolic equations with piecewise contin- 
uous initial conditions is developed. The new scheme is proven to be stable in the energy 
norm for both continuous and discontinuous solutions. In contrast to the existing high- 
resolution shock-capturing schemes, no assumption that the reconstruction should be total 
variation bounded (TVB) is explicitly required to prove stability of the new scheme. A 
rigorous truncation error analysis is presented showing that the accuracy of the 3rd-order 
ESWENO scheme is drastically improved if the tuning parameters of the weight functions 
satisfy certain criteria. Numerical results show that the new ESWENO scheme is stable 
and significantly outperforms the conventional third-order WENO finite difference scheme 
of Jiang and Shu in terms of accuracy, while providing essentially nonoscillatory solutions 
near strong discontinuities. 


I. Introduction 

Higli-order finite difference weighted essentially nonoscillatory (WENO) schemes have now become a 
state-of-the-art methodology for solving problems with strong discontinuities that arise in computational 
fluid dynamics, computational aeroacoustics, computational electromagnetics, and other areas. The WENO 
schemes have evolved from the successful ENO reconstruction idea which was originally proposed by Harten. 1 * 
In contrast to the ENO reconstruction which adaptively chooses one ’’smoothest” stencil from p candidate 
stencils, the WENO schemes use a weighted combination of all the candidate stencils to form the WENO 
reconstruction of the flux. Each weight is a nonlinear function of so-called smoothness indicators that measure 
the smoothness of the numerical solution on that particular stencil which this weight has been assigned. If 
the solution is smooth in all candidate stencils, then the WENO scheme asymptotes to a target scheme 
which provides either higher order accuracy 2,3 or improved resolution for short waves. 4,5 If the solution 
is discontinuous in one of the candidate stencils, the WENO scheme effectively nullifies the corresponding 
weight, thus biasing the stencil away from the discontinuity and successfully emulating the ENO stencil 
choosing procedure. To date, various WENO schemes have been developed for both finite difference, 3-7 and 
finite volume 2, 8-10 formulations. 

Although the higli-order WENO schemes have been successfully used for solving a broad spectrum of 
applied problems, very few theoretical stability proofs for these schemes are available in the literature. It is 
shown in 3 that the WENO scheme of Jiang and Shu is stable for the conservation law equation u t + V-f(u) = 
0, if a smooth flux splitting and a nth-order Runge-Kutta scheme (n > 3) are used, and the differential 
problem is at least p+ [(d + 1) / 2] + 2 times continuously differentiable (d is the dimensionality of the problem, 
and p is the order of approximation). For the scalar conservation law equation in one dimension, Jiang and 
Yu 11 have shown the existence and stability of discrete stationary shocks for the 3rd-order WENO scheme in, 3 
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based on the global Lax-Friedrichs flux splitting. Note, however, that for the 3rd-order WENO scheme, the 
key condition for existence and stability (see Theorem 4.4 in 11 ) has been verified numerically only for sample 
values of the parameter qo, rather than proven rigorously for all qo € [0, 1]. Recently, Ferretti 12 as well as 
Qiu and Shu 13 proved the convergence of a quite general class of self-similar total variation bounded (TVB) 
schemes to the entropy solution of strictly convex scalar one-dimensional Hamilton- Jacobi and conservation 
law equations under the large time step assumption: At/ Ax — > +oo. This assumption allows one to 
eliminate the need for the cell 14 or wavewise 15 entropy inequality that is traditionally required to prove the 
convergence of this class of TVB schemes. Note that this time step constraint violates the CFL condition, 
a natural condition for time-dependent hyperbolic problems (e.g. fluids, acoustics, electromagnetics). Thus, 
these TVB schemes are only amenable for implicit temporal integrators. 

In contrast to nonlinear ENO and WENO formulations, stability proofs for linear schemes are elementary. 
Indeed, Kriess and Scherer 16 recognized that discrete operators satisfying a summation-by-parts (SBP) 
condition, are automatically stable in an L 2 -energy norm. Furthermore, they noted that the SBP property 
could be used to develop new schemes, not just prove the stability of existing ones. Numerous schemes 
currently exist that satisfy the SBP property. In addition to weak form FEM operators, SBP operators exist 
for central schemes, 1 ' Pade schemes, 18 upwind schemes, 19 and spectral collocation schemes. 20-22 Although 
proofs of stability exist for linear schemes, robustness is frequently lacking. Because they use a fixed stencil, 
they are susceptible to Gibbs oscillations for discontinuous data. 

Drawing from these two complementary approaches (WENO and SBP), we develop a new third-order 
WENO-type scheme (called ESWENO) suitable for linear wave equations with discontinuous solutions, and 
prove that the new scheme is Energy Stable, i.e. stable in an E 2 -energy norm. By construction, the ESWENO 
scheme satisfies nonlinear SBP and positive semidefiniteness conditions at each instant in time, which provide 
stability. Thus, E 2 strict stability is attained without the need for a TVB flux reconstruction, or a large-time- 
step constraint. Despite the SBP stencil constraint, the new ESWENO scheme still retains the underlying 
WENO character of the background scheme. Indeed, numerical experiments demonstrate that the ESWENO 
scheme delivers stable essentially nonoscillatory solutions for problems with strong discontinuities. 

This paper is organized as follows. In Section 2, an energy estimate for the continuous singular perturbed 
wave equation is obtained, followed by an energy estimate for the corresponding discrete problem. In Section 
3, we show that the dissipation operator of the conventional 3rd-order WENO scheme is not positive semidef- 
inite, thus indicating that the scheme may become locally unstable if strong discontinuities or unresolved 
features are present in the domain. In Section 4, we present the 3rd-order ESWENO scheme and prove that 
the new scheme is stable in the energy norm, 3rd-order accurate for smooth solutions including local smooth 
extrema, and conservative. Furthermore, based on our rigorous truncation error analysis, we show how to 
choose the tuning parameters of the weight functions to significantly improve the dissipation properties of the 
ESWENO scheme. In Section 5, we extend the scalar ESWENO scheme to hyperbolic systems and obtain 
an energy estimate for the characteristic form of the hyperbolic system of equations, while a similar estimate 
is not available if the system is discretized in the component-wise fashion. In Section 6, we present numerical 
experiments that corroborate our theoretical results. We summarize and draw conclusions in Section 7. 

II. Energy Estimates 

Semi-discrete summation- by-parts (SBP) operators are constructed to explicitly satisfy an E 2 -energy 
estimate. That is, the semi-discrete operator ’’mimics” the continuous operator in terms of the energy of the 
system. We begin by presenting a derivation of the energy estimate for the continuous system, followed by 
a derivation of the mimetic SBP semi-discrete operator. 

II. A. The Continuous Problem 

Let us consider a linear, scalar wave equation with periodic boundary conditions 

§y + §7 = °> / = au > t>0, 0<x<l, 

u(0, x) = Uo(x) (1) 

u(t, 0) = it(f, 1) 
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where a is a constant, and u o(x) is a bounded piecewise continuous function. Without loss of generality, we 
assume that a > 0. 

Along with Eq.(l), we also consider the following singular perturbed wave equation which is subject to 
periodic boundary conditions 


f + i = / = ««> *> 0 , 0<*<1 

u( 0, x) = uq(x) 


(2) 


where B is a nonlinear positive semidefinite differential operator that depends only on u and its derivatives, 
i.e. , (v,Bv) = Jq vBvdx > 0 for all sufhciently smooth functions v{x) that satisfy the boundary conditions. 
The right hand side of Eq. (2) can be interpreted as an artificial dissipation term for a finite difference 
scheme approximating Eq. (1). 

To obtain the energy estimate, Eq. (2) is multiplied by u and integrated over the entire domain, which 
yields 




uBu 


■®io 


l 

J u x Bu x dx, 
o 


(3) 


where || • ||; 2 is the continuous 1-2 norm. By taking into account that the boundary conditions are periodic 
and B depends only on u and its derivatives, the second term on the left hand side and the first term on the 
right hand side of Eq. (3) vanish, and the following energy estimate can be obtained 


1 

d / • 

^IMI? 2 = -2 / u x Bu x dx < 0 (4) 

0 

Because B is a positive semidefinite differential operator, the right hand side of Eq. (4) is nonpositive, 
thus providing the required dissipation property. It should be emphasized that in the above derivation, 
no additional assumptions have been made about B other than its positive senridefiniteness and that the 
operator B depends only on u and its derivatives, which makes Eq. (2) nonlinear. 


II. B. The Discrete Problem 

The energy estimate derivation for equation (2) relies on two properties of the continuous operator. First, 
the derivative operator must satisfy the integration-by-parts (IBP) property, i.e., (v,f x ) = — (v x ,f) (for 
the periodic case). Second, the function B must be positive semidefinite, i.e. (v,Bv) = vBvdx > 0. In 
the present analysis, a stable, semi-discrete, third-order finite difference scheme for Eq. (1) is constructed by 
approximating f x as a sum of two terms. The first term mimics the IBP property, while the second mimics 
the senridefiniteness of the operator B. This scheme is presented next. 

Define a uniform grid x 3 = j Ax, j = 0, J, with Aai = 1 / J. On that grid, define a flux f = au and its 
derivative f x = au x , where u = [u(xo,t), . . . ,u(xj,t)] T and Uj, = [u x (xo,t), . . . , u x (xj, t)] T are projections 
of the continuous solution and its derivative onto the computational grid. Finally, define an approximation 
for the first-order derivative term in Eq. (1) as follows: 

— = (D c + D a )i + 0(Ax 2p ~ 1 ), (5) 

The matrix D c is a nonlinear central finite difference operator given by 


D C = P~ 1 Q ; f x - F-iQf = 0(Ax 2 P) 
P = AxI ; Q + Q T = 0, 
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where I is a ( J + 1) x ( J + 1) identity matrix. Note that the approximation (6) satisfies the SBP rule 

(u, D c v) P = ~(D C u, v) P , (7) 

with (u,v)p = u T P\. Thus, the first mimetic relationship between the continuous and discrete equations 
is satisfied by D c \ i.e. IBP — > SBP. (The relation P = A xl is valid only for the periodic case. Alterations 
in P result for the nonperiodic case, but this generality is outside the scope of the current work.) 

The matrix D a is a nonlinear finite difference artificial dissipation operator, which mimics its continuous 
counterpart ^ (P^). It is given by 


D a = p-^Df SD 1 ; p- 1 DjSD 1 u = 0(Ax 2p - 1 ) ; v T {S + S T )v> 0 (8) 

where P is the same positive definite matrix used for the approximation of the first-order derivative term Eq. 
(6), v is an arbitrary real- valued vector of length ( J + 1), and D\ and D f are two-point difference operators 
given by 


II 

Q 

-1 1 

o 

II 

1 

1 

o 

1 — 1 
1 


o 



o 

. J 


Equation (8) requires that the matrix S must be positive semidefinite, and that the nonlinear artificial 
dissipation term should be zero to the design order of the scheme (at least in regions where the solution of 
the differential problem (1) is sufficiently smooth). 

Using the SBP operators Eqs. (6) and (8), Eq. (1) is discretized as follows: 

— + P-'Qf = -P" 1 ^ SD 1 f, (9) 

where f = au and u = [uo(t) , Ui(t) , . . . ,uj(t),] T is the discrete approximation of the solution u of Eq. (1), 
Q and S are nonlinear matrices, i.e., Q = <3(u) and S = S'(u). To show that the above finite difference 
scheme is stable, the energy method is used. Multiplying Eq. (9) with u T P yields 

5dp U llp + ou t Qu = -a{Diu) T SDiu, (10) 

where || • ||p is the P norm, i.e., ||u||p = u T Pu. Adding Eq. (10) to its transpose, we have 

|||u|| 2 P + au T (Q + Q t ) u = -a (Pmf (S + S T ) D lU . (11) 

Taking into account that the boundary conditions are periodic and Q is fully skew-symmetric, the second 
term on the left hand side vanishes, and the energy estimate becomes 

^||u||! = -a (P lU ) T (S + S T ) D lU < 0. (12) 

The right hand side of Eq. (12) is nonpositive, because the matrix S is positive semidefinite and a > 0, thus 
providing stability of the finite difference scheme Eq. (9). This result can be summarized in the following 
theorem: 

Theorem 1. The approximation (9, 6, 8) of the problem (1) is stable if Eqs. (6, 8) hold. 

Remark 1. Despite that the initial boundary value problem (1) is linear, the finite difference scheme (9) 
constructed for approximation of Eq. (1) is nonlinear, because it is assumed that the matrices Q and S (and 
in principle P) depend on the discrete solution u. 
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Remark 2. The only constraints imposed on the matrices Q and S are skew-symmetry of the former and 
positive semidefiniteness of the later. No other assumptions have been made about a specific form of the 
matrices Q and S to guarantee the stability of the finite difference scheme (9). 

Remark 3. The discrete operators defined by equations (5 - 8) are symbolically identical to those used for 
conventional SBP operators (see 15, 16, 22 ). What is new, however, is the fact that the matrices Q and S 
depend on u. 


III. Conventional 3rd-order WENO Scheme 


III. A. Definitions 

The conventional WENO scheme 3 for the scalar 1-D wave equation (1) with a > 0 is given by 

duj 


fW _ fW 
•0+5 0-5 


dt. 


Ax 


= 0 


(13) 


where L is the 3rd-order WENO flux which is 
2 


(15) 


fW — / ,0 _i_ / .1 fl 1 ) (AA\ 

Jj+i ~ U j+l/2Jj+l/2 + ^.7+1/2 «C?‘+l/2 ’ ( 14 ) 

where u>° and uj 1 are weight functions assigned to two stencils {xj,Xj+i} and {xj-i,Xj}, respectively. The 
2nd-order fluxes, r = 0, 1, in Eq. (14) are defined on these two different stencils as follows: 

fj+ 1/2 = + hf(uj+ 1), 

fj+ 1/2 = 1 ) + |/(«j) 

In the present analysis, if the flux in Eq. (1) does not possess the property f'{u) = a > 0, then the global 
Lax-Friedrichs flirx splitting is used 

/ ± (m) = ^ (/(«) ± A max u) , (16) 

where A max is a constant that satisfies the following constraint 

Amax ^ I & I. 


(17) 

with respect to a point j + \- 
The weights and smoothness indicators used herein are very similar to those proposed by Jiang and Shu 
in, 3 which are given by 


The 2nd-order stencils for /\_ 1 are a mirror image of the stencils used for 


w . , i = — -f — , 

] + ^ “0+«l ’ 


cy,f — 
OLy — 


F+^r)” 

d r 


r = 0,1 


with 


j.q 35 do 3 ? 0o On 4 1 Uj ) 


do = 


di = 


di = 


01 = (Uj - Uj-!) 


(18) 


(19) 


where k and m are positive constants, e is a small positive parameter that does not allow the denominator 
to become equal to 0. Note that if k = 2, m = 2, and e = 10~ 5 6 , then the weights and smoothness indicators 
given by Eqs. (18, 19) are identical to those used in. 3 The weights (18, 19) possess the following two 
properties 


0 < w] +1 < 1, 




■ w- 


= 1, 


W r . 1 = dr 
d“T 2 

w r 1 = d r 
J-5 


3+ 

0( Ax k+1 ) 
0(Ax k+1 ) 


for r = 0,1, j = 0, J 
for r = 0, 1 


(20) 

(21) 


Note that Eq. (20) is valid for any vector u, whereas Eq. (21) holds only where the solution is sufficiently 
smooth. 
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III.B. Stability Analysis 

Before constructing a new energy stable WENO (ESWENO) scheme, we first show that the conventional 
3rd-order WENO scheme developed in 3 can be represented in the form of Eq. (9) with appropriate choices 
for the matrices P, Q , and S. Note, however, that the matrix S of the WENO scheme does not, satisfy 
all the sufficient conditions given in Eq. 8. Thus, it is not immediately evident how to form an the energy 
estimate for this scheme. 

To show this, we begin by defining the 3rd-order WENO derivative matrix 


Pwf 


f w 

o+i 


fW 
Jj-i 
j 2 


Ax 


(22) 


Substituting Eq. (14) and (15) into (22) (and without loss of generality assuming f'(u) = a > 0) yields the 
periodic pentadiagonal matrix Dw, 


D^y — -P 67 it(i [dj — 2 5 dj — r , dj , djy \ , dy-j-2] 


with entries 


dj- 2 

T 

W 1 1 
3-2 

dj- 1 


-3 W )_,-w) H -w«_ h 

dj 

1 

2 Ax 

3 w j+i +w° j+ i 

dj+ 1 


W° 1 

j 


3+i 

L «i +2 J 


0 


(23) 


If the symmetric part of Dw is positive semidehnite for all values of the parameters w * (thereby emulating 
the conditions of Theorem 1) then stability immediately follows. To represent Dw in the form of Eq. (9), 
we decompose Dw into its symmetric and skew-symmetric components 


% = 2 I-^w ~ Pwf] + 2^ w + = D w &W + D w’ n - 


(24) 


Equating the skew-symmetric and symmetric components of Dw with the equivalent terms from Eq. (9) 
yields 


D sk ew = P iq = Penta[qj-2,qj-i,qj,qj+i,qj+2] 

D syrn _ p-i^T^w^ = Penta [rj - 2 , rj-i , rj, Tj+1 , r j+2 ] 


(25) 


with 


and 


_ _ 

T 

* 

Qj- 2 


•? 2 



— 3ud 1 — w 1 . 1 — 2w° 1 

Qj- 1 


3 ~ 5 J+5 J-5 

Qj 

1 

4Aai 

0 

Qj+ 1 


3^ +| +^ + 3+2 % ° +| 

- Qj + 2 - 


— tu 1 , 3 



L J+f J 


’ r j-2 ’ 
V 3 - 1 

T 

" 

j 2 

r j 

1 

AAx 

2(3^ 1 +| -^°_ | +^° +| ) 

r i+ 1 



- T?+2 . 


1 

A 

+ 

to|co 

1 


-1 T 


(26) 


( 27 ) 
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The pentadiagonal matrix Q does not contribute to the P-norm of the energy, because it is entirely 
skew-synnnetric. Nevertheless, inspection of Eq. (26) immediately reveals that the matrices P and Q are 
given by 

P = A xl ; Q = PD$ ew . (28) 

The diagonal P matrix is symmetric positive definite, making it suitable as a norm. 

Let us show that the symmetric part of Dw is not positive semidefinite. Lengthy manipulations of the 
expression 

D svm = p-i]jT S w Di _ p en t a [ rj _ 2 , rj _ 1 , rj , rj+1 , rj+2 ] 

reveal that periodic matrix S w derived from the conventional 3rd-order WENO scheme (13-15, 18,19) is 
tridiagonal and given by the form 


S w = Tri 



1 

3 - 


1 > 

2 


2 W J-5’ 




(29) 


Given Eq. (29), it is still not apparent whether or not the symmetric part of S w is positive semidefinite, 
the sufficient condition for stability of the discretization (9). Further tedious manipulations reveal that the 
tridiagonal matrix S w in the WENO artificial dissipation operator can be written as 


S w = S? +D‘?S% f D u 

where S ^ and S ^ are the following diagonal matrices: 



(30) 


(31) 


From Eq. (31) it immediately follows that Sw in Eq. (29) is not positive semidefinite and therefore 
does not satisfy the conditions of Theorem 1. Indeed, the diagonal matrix S™ is positive semidefinite 
because w\ i > 0 for all j, while the diagonal matrix S^ given by Eq. (31) is not positive semidefinite 
because s™.. can be either positive or negative. As a result, the 3rd-order WENO scheme of Jiang and Shu 
does not provide the energy estimate obtained in Section 2. Furthermore, near strong discontinuities or 
unresolved features, s™ may become 0(1), thus shifting the eigenvalues of the WENO artificial dissipation 
operator —P~ 1 DjS w D i to the right half plane. Because Theorem 1 provides only the sufficient condition 
for stability, no conclusion can be made about stability of the conventional WENO scheme for Eq. (1) with 
a discontinuous initial condition. Note that in, 3 the stability of the WENO scheme has been proven only for 
sufficiently smooth solutions. 


IV. 3rd-order Energy Stable WENO Scheme for the Scalar Wave Equation 

We now construct the operators P, Q, and S in Eq. (9) such that all the conditions of Theorem 1 are met. 
As has been shown in the foregoing section, the operators P and Q defined by Eq. (28) satisfy the conditions 
of Theorem 1 and can be directly incorporated into the new ESWENO scheme. The only condition that 
has not been met is positive semidefiniteness of the matrix S in the artificial dissipation operator D a . We 
propose the following symmetric matrix S in Eq. (9), which possesses the required properties (8): 

S = S 1 +D'[S2D 1, (32) 

where S± and S 2 are the following diagonal matrices: 

s hd = /h 

^ =H-i 

where 

N = \^ w ) +h -w)_ h ) 2 +P, 


(33) 

(34) 
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and 6 is a small positive parameter which has to be chosen so that the ESWENO artificial dissipation term 
is 3rd-order accurate. Possible ways how one can chose the parameter S are discussed in Section 4.1. Note 
that fij is a C°° function of u, which is consistent with the smoothness of the weight functions. 

The ESWENO scheme is given by Eq. (9) with P, Q, and S defined by Eqs. (28, 32-34) and the 
weights (18-19). We will show next that the ESWENO scheme is: 1) stable in the energy norm, 2) 3rd-order 
accurate, if the solution of Eq. (1) is sufficiently smooth and some constraints are imposed on k and <5, and 
3) conservative. 

Remark 4. The 3rd-order ESWENO scheme has a 5-point stencil which is as wide as the stencil of the 
conventional 3rd-order WENO scheme. 

IV. A. Stability 

Theorem 1 provides the sufficient conditions (6, 8) for stability of the finite difference ESWENO scheme 
(9). Therefore, if the constraints (6) and (8) are met, then the proposed ESWENO scheme is stable in the 
energy norm. It can easily be seen that the matrices P and Q defined by Eqs. (28) satisfy the conditions (6). 
Therefore, to prove the stability of the proposed ESWENO scheme, one only has to show that the real- valued, 
nonlinear operator S given by Eqs. (32-34, 18, 19) is positive semidefinite, i.e., v T S'v = v T ( S+ 2 S ) v > 0 

for all discrete real- valued vectors v of length (J + 1). Let us show that the matrix S constructed in the 
foregoing section is positive definite, thus satisfying the conditions of Theorem 1. Indeed, the diagonal matrix 
S 2 given by Eq. (33) is positive semidefinite, because lu 1 _ 1 > 0 for all j. From the following inequality 

which holds for any weight function w 1 i if S > 0: 

3~r 2 

it immediately follows that each diagonal element si.. is positive, thus providing positive definiteness of Si. 
Hence, the symmetric matrix S given by Eq. (32) is positive definite. Note, however, that SD\ is positive 
semidefinite rather than positive definite, because for v = (1, . . . , 1) T , we have (Div) T SDiV = 0. It should 
also be noted that the matrix S depends on the discrete solution u of Eq. (9) and is positive definite for any 
real-valued vector u regardless of whether u is continuous or discontinuous. 

Remark 5. The parameter fij can be chosen in many different ways so that it provides the energy estimate 
and consistency. Among possible candidates for /j.y , there is at least one that requires no tuning parameters, 
which is given by 

to = 4 ■ 

The main drawback of ji 3 presented above is that the additional “dissipation” term P _1 Pf diag(/zi , . . . , 
which provides the energy estimate, is not actually dissipative and therefore makes the ESWENO scheme 
more susceptible to spurious oscillations. 

Remark 6. Only the property (20) has been used to prove the positive definiteness of the matrix S. This 
property does not depend on a particular form of the weight functions and remains valid for any weight 
functions that satisfy Eq. (20). 

Remark 7. The energy estimate cannot be obtained if the weights are negative, because in this case, the 
matrix S is not positive semidefinite. 

IV. B. Consistency 

We will now show that for sufficiently smooth solutions, the ESWENO scheme given by Eq. (9, 28, 32-34) 
with the weight functions (18, 19) is third-order accurate if the parameters k in Eq. (19) and <5 in Eq. (34) 
satisfy the following constraints: 

k > 0 (35) 

<5 = 0(J~ 2 ), (36) 
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where J is the number of grid cells. Note that the above constraint for S does not depend on the size of the 
computational domain, thus providing scale invariance of the ESWENO scheme. In other words, the scheme 
is invariant when the spatial and time variables are scaled by the same factor. 

In, 3 it was proven that the conventional WENO scheme (13-15) is 3rd-order accurate if the weights 
(18-19) satisfy the following constraints: 


tu r = d r + 0( Ax) 
u r = d r + 0( Ax), 


(37) 


and the solution is sufficiently smooth over all candidate stencils. As will be proven next, if k in Eq. (19) is 
strictly positive and the solution of the continuous problem is sufficiently smooth, then the weight functions 
satisfy Eq. (37). Indeed, expanding the smoothness indicators Po and Pi defined on stencils {xj+i+i , Xj+i } 
and {xj+t-i, Xj+i} in the Taylor series about Xj, we have 

Po = (uj+i+i - Uj+i) k = (u x ) k Ax k + fc(2 2 +1) (u x ) k ~ 1 u xx Ax k+1 + 0{ Ax k+2 ) , . 

Pi = ( Uj+i - Uj+i- if = ( u x ) k Ax k + fc(2 2 ~ 1} ( u x ) k ~ 1 u xx Ax k+1 + 0( Ax k+2 ), 

where l is an arbitrary integer number, Uj = u(xj), u x = u x {xj ), and u xx = u xx {xj) are projections of the 
continuous solution and its derivatives on the computational grid. Substituting Eq. (38) into Eq. (18) yields 


U) 


j+(+§ 


do~\-di 


do+di 


do 

’ l + (flx) fc ^A+ M2 j + 1 ) (Cx) fc - 1 flxxAxfc + l+o(Ax|j±) 

l + (5 I )t4fl4.M^zl)(5,)t- 1 tl ,xAxHl + o( AJ') +2 ) 

do 

(l+^(«x) fc - 1 MxxAx'=+i+o(^±^)) m 


(39) 


Equation (39) has been derived under the following assumption: 

Ax k < e, (40) 

which is valid on sufficiently fine grids. Expanding the denominator of Eq. (39) about the same point Xj 
and using do = 2/3 and d\ = 1/3, we have 


w° +;+ i = d 0 - 


2kmAx k+1 

9e 


u x ) 


k- 1 = 


o 


^ Ax k + 2 


(41) 


Note that the leading truncation error term of to® +l+ i at Xj does not depend on l , which implies that it is 
invariant to an arbitrary shift of the stencil. The leading truncation error terms of the other weight functions 
can be obtained by using the same procedure outlined above, thus leading to 


“l+l+h = 


{ _ l)r 2kmA^ ^-1^ + 0 {A^A) 
(-l) r 2fcm t xfc+1 + O(^), r = 0/T. 


(42) 


As follows from Eq. (42), if the parameter k is strictly positive then the conditions (37) are met, which 
indicates that the conventional WENO scheme with the modified weight functions defined by Eqs. (19- 
20) is 3rd-order accurate in regions where all required derivatives of the solution are continuous. Another 
important conclusion that can be drawn from Eq. (42) is that the leading truncation error terms of u> r and 
Co r , r = 0, 1 are equal to each other and independent of the stencil on which they are defined. Note that k in 
Eq. (19) cannot be equal to 0, because in this case, the weight functions given by Eqs. (19-20) degenerate 
to constants, thus leading to the 3rd-order, linear upwind-biased scheme that does not possess the ENO 
property. As will be shown in the next section, a stronger constraint should be imposed on the parameter k 
to guarantee that the ESWENO scheme is 3rd-order accurate near smooth extrema. 

The P and Q operators of the ESWENO scheme are identical to those of the conventional 3rd-order 
WENO scheme. Comparing Eqs. (30-31) and (32-34), one can see that the ESWENO scheme can be 
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obtained from the WENO scheme of Jiang and Shu by adding the following additional artificial dissipation 
term: 


[-Desii - flwu] • = 


- Mj+i 


V . , Q —UJ . 


( fj - fj+ 1)» 


(43) 


where /i ? is defined by Eq. (34) and Z?es and Dw are the ESWENO and WENO operators, respectively. 
Let us show that for sufficiently smooth solutions, the additional ESWENO dissipation term (43) is of the 
order of 0( Ax 3 ) if Eqs. (35, 36) hold. Expanding Eq. (43) in the Taylor series and taking into account Eq. 
(42) yield 


Ax 


4 + o (*$£■) + o (^i 
0(5Acc) + O ( A ^- +J ) + O ( 


Ax 


4Ax 

k + 3 N 


(44) 


The above equation has been derived using fj + 1 — 2 fj + fj - 1 = 0(Ax 2 ). From Eq. (44) it immediately 
follows that the parameter <5 has to satisfy the criterion (36) to guarantee that the ESWENO scheme is 
third-order accurate for sufficiently smooth solutions. As a result, Eq. (44) can be recast as follows: 


^'+2 


I — e 


3- ; 


Ax 


= O 


(A x 2k+3 \ 






+ o 


Ax k+3 ^ 


(45) 


For k > 0, the right hand side of Eq. (45) is of the order of 0(Ax 3 ), thus providing the design order of 
accuracy for the ESWENO scheme. Note that for the weights (18, 19), the order of approximation of the 
ESWENO scheme does not depend on the parameter m in Eq. (18). 


IV. C. How to Choose Parameters k, m, and e 

In contrast to the original WENO scheme 3 which has only one tuning parameter e, the ESWENO scheme 
with the modified weight functions (18, 19) and the additional artificial dissipation term (43, 34) has four 
tuning parameters: k. m, e, and S. In connection with this, a question that needs to be addressed is how to 
choose these parameters to improve the accuracy of the ESWENO scheme while preserving the nonoscillatory 
properties of the original WENO scheme near strong discontinuities. As has been shown in the previous 
section, the parameter S should satisfy Eq. (36) to provide the design order of accuracy, thus leaving only 
three tuning parameters that have to be determined. 

As will be shown in Section 6, the conventional 3rd-order WENO scheme is too dissipative as compared 
with the corresponding 3rd-order, linear, upwind-biased scheme which is obtained by setting u> r = d r and 
uj r = d r , r=0, 1. The main reason for such a behavior is the following choice of the tuning parameters: 
k = 2, m = 2, and e = 10~ 6 used in the conventional WENO scheme. Furthermore, for this choice of the 
parameters, the conventional 3rd-order WENO scheme may locally exhibit only a 2nd-order convergence 
rate even on moderate grids that are widely used in practical applications. Indeed, using Eq. (38) with 
k = 2, m = 2 and assuming that e = 0(Aa: 3 ) on a given mesh, uA +i at Xj can be evaluated as 


do~hdi 

e+(u x Ax+±u xx Ax 2 +0(Ax 3 )^ 2 


_ e +(oArr-^O xx Ax2 + 0(Ax3)) 2 _ 



do 


do+di 

[ e/Ax 2 +u 2 +u x u xx Ax + 0(Ax 2 ) 2 


[ e/ Ax 2 + u 2 — u x u xx Ax + 0(Ax 2 ) 



(46) 


For a given Ax, one can come up with an example when there exists a point x* such that either e/Aar +u 2 + 
u x u xx Ax = 0 or e/Aar +it 2 — u x u xx Ax = 0 near local extrema of the continuous solution, where u x decreases 
and passes through 0 and u xx = 0(1). Considering the above relations as quadratic equations for u x , it can be 
easily shown that these equations have real- valued solutions if \u xx \ > , thus indicating that the point x* 
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exists for sufficiently coarse grids. Without loss of generality, one can assume that e/ Ax 2 + u 2 +u x u xx Ax = 0 
at Xj = x*. As a result, Eq. (46) can be recast in the following form: 


0+2 


do~\-d\ 


do 

Q( Ax 2 ) 

2e/ Ax' 2 +2u'£+0(Ax' 2 ) 


do 

do+0(Ax 2 ] 


= 1 + 0( Ax 2 ) 


(47) 


Equation (47) has been derived under the assumption that e = 0( Ax 3 ) which can be interpreted as a 
constraint on the grid spacing Ax for which Eq. (47) holds. From Eq. (47), it follows that w° +1 / 2 ^ do and 
the conventional WENO scheme may locally exhibit only the 2nd-order convergence rate if the grid spacing 
is of the order of 0(e 1//3 ). For e = 10 -6 , Aa’ = 0{e 1 ^) = O(10 -2 ) which indicates that this deterioration in 
accuracy may occur even on moderate grids. 

Remark 8. Estimate (47) holds only for sufficiently coarse meshes: Aa: > 0(e - 1 / 3 ), whereas asymptotically, 
the conventional WENO scheme is third-order accurate near smooth extrema. 

With the same tuning parameters used for the conventional WENO scheme, the ESWENO scheme 
becomes even more dissipative than its counterpart, because of the presence of the additional artificial 
dissipation term (43). The accuracy of the ESWENO scheme (and WENO) can be drastically improved if 
the parameters k, m, and e are chosen based on criteria discussed next. 

The truncation error analysis Eqs. (38-44) suggests that the grid spacing Ax in Eq. (18) should satisfy 
the following constraint: e ^ Ax k . As follows from Eq. (45), the contribution that comes to the leading 

truncation error term from the weight functions is of O 1- Ax e +3 ^ O (Aa; 3 ). Hence, if the condi- 

tion (40) is satisfied, then the ESWENO scheme can provide a 3rd-order convergence rate on meshes with 
Ax = 0(e 1 ^ k ), which can be controlled by an appropriate choice of the parameter e. Equation (40) gives an 
estimate for the largest grid spacing on which one should expect to obtain the design order of convergence 
of the ESWENO scheme. It should also be emphasized that Eq. (40) bounds the grid spacing from above, 
so that if this criterion is met on a given grid, it will also be satisfied on any finer mesh. 

Another important observation is that the ESWENO scheme remains 3rd-order accurate even near smooth 
extrema if Eq. (40) holds and k > 1. Indeed, using Eqs. (42) and (40), w° A near local extrema can be 
evaluated as follows: 

0 2kmAx k+1 xfc_i_ _.,Aa; fc+2 

1 = d r ( u x ) k 1 u xx + 0( ) (48) 

2 9e g 

At the smooth extrema, u xx is finite, while (■u a; ) fc_1 becomes infinitely large if k < 1. Therefore, to keep the 
(u x ) fc_1 term bounded near the smooth extrema, the following constraint should be imposed 


k > 1. 


(49) 


If the conditions (40) and (49) are met, Eq. (48) becomes 

i = d r + 0( Ax), (50) 

which is sufficient to guarantee that the WENO scheme is third-order accurate, as has been shown in. 3 
Moreover, the additional ESWENO dissipation term also remains third-order accurate near smooth extrema 
if k > 1 and e Ax k , as follows from Eq. (45). Hence, if Eqs. (40, 49) hold, then the ESWENO scheme is 
3rd-order accurate in regions where the continuous solution is sufficiently smooth, including smooth extrema. 

A question that has not been addressed yet is what constraints should be imposed on the tuning pa- 
rameters e, k, and m to emulate the ENO stencil biasing strategy. Let u{x) be a piecewise smooth function 
on an interval [xj- 2 , %j+ 2 \ such that the discontinuity is located at Xd : Xj < Xd < Xj+i- Then, the weight 
function u>® +1 / 2 can be evaluated as follows: 


A* 


do + d\ 


do 

Q(i) 

e+0( Ax k ) 


O(l) (e m + e m-1 0(Ax fc )) 


(51) 
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To emulate the ENO property near the strong discontinuity, the stencil Sj = {xj,Xj+ 1 } has to be eliminated 
from the approximation, which can be achieved if the corresponding weight satisfies the following constraint: 


w j+1/2 ^ 1 - 

Taking into account Eqs. (40) and (51), the constraint (52) becomes 

^+i/ 2 = 0(e m ) « 1. 


(52) 

(53) 


Note that the smaller the weight, w® +1 ^ 2 , the stronger the stencil is biased away from the discontinuity, thus 
locally reducing the order of approximation and introducing more dissipation. To satisfy Eq. (53) , it suffices 
to have 


f«l 

m > 0 


(54) 


Summarizing, we can conclude that if the constraints (36, 40, 49, 54) are met, then the ESWENO 
scheme is 3rd-order accurate in regions where the solution is smooth and provides essentially nonoscillatory 
solutions for problems with strong discontinuities. Another observation is that the weights should be at least 
C 4 functions of the solution to guarantee that the ESWENO scheme is 3rd-order accurate. A good choice 
for the parameter k is k = 2, which both satisfies Eq. (49) and provides that the weights given by Eqs. (18, 
19) are C°° functions of u. Note that it is not desirable to use k larger than 2 because it makes the scheme 
more prone to spurious oscillations. 

The above consideration also suggests that e should be chosen much greater than 10 “ 6 used in the 
conventional WENO scheme, so that the constraints (40) and (54) hold. Equation (40) gives an estimate for 
the maximum Ax for which the ESWENO scheme should exhibit the design order of convergence. Based 
on our numerical experiments presented in Section 6, a reasonable choice for e is 0(1O -3 ). This choice of e 
and k translates into the following estimate for Ax : Ax < 0( 10~ 3//2 ) , thus indicating that the design order 
of accuracy can be achieved on coarse grids. Another way of choosing the parameter e is to represent it as 
a function of the number of grid points J, so that the constraints (40) and (54) are met. A possible choice 
for e is 

e = 0{J~ l ). (55) 

The parameter e given by Eq. (55) depends only on the number of grid points and is independent of the size 
of the computational domain, thus preserving the scale invariance of the ESWENO scheme. For this choice 
of e and k = 2, the constraint (40) reduces to 

0(Ax 2 ) « 0(J- 3 ), 


which is valid if the number of grid points is much greater than L 2 , where L is the nondimensional length of 
the physical domain. Therefore, Eq. (55) can provide the design order of accuracy of the ESWENO scheme 
even on coarse grids. 

Our choice of the parameter m is determined by the balance between the accuracy in smooth regions 
and good shock-capturing capabilities. On the one hand, Eq. (42) shows that the leading truncation error 
term is directly proportional to m. Therefore, a reduction of the parameter in decreases the coefficient in 
front of the leading truncation error term, thus making the ESWENO scheme less dissipative. On the other 
hand, the parameter m cannot be set too close to zero, because in this case, Eq. (53) imposes a very severe 
constraint on the parameter e : e m <C 1. In all our numerical computations, we use rn = 1, which provides 
good dissipation properties both in smooth solution regions and near strong discontinuities. Note that the 
above conclusions are valid not only for the ESWENO scheme, but also for the conventional 3rd-order WENO 
scheme. 


IV. D. Conservativeness 

ENO and WENO schemes have emerged as methods that provide both high-order accuracy and excellent 
shock-capturing capabilities. These properties are extremely important for solving problems that are char- 
acterized by the presence of a lot of structure in the smooth part of the solution, strong discontinuities, and 
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their interaction. As follows from the Lax-Wendroff theorem, 24 if a conservative numerical scheme approx- 
imates Eq. (1) and converges, then it converges to a weak solution of the conservation law equation (1). 
Therefore, to accurately calculate the strength and speed of shock waves and contact discontinuities with 
finite mesh size, a numerical scheme should be conservative at the discrete level. 

For the continuous problem Eq. (1), integration over the spatial domain leads to 

l 

J udx + f(t , 1) - f(t, 0) = 0. (56) 

o 


The above equation shows that the total quantity of a conserved variable u in the domain changes only 
because of the flux through the domain boundaries. As has been shown in Section 4.2, the ESWENO 
scheme can be written in the semi-discrete conservative form (13) with the flux / ES = / w + e, where / w is 
the conventional WENO flux, and e is the additional dissipation flux term given by Eq. (43). From these 
equations it immediately follows that the numerical flux telescopes across the domain to the boundaries, thus 
mimicking the conservation property of the continuous problem and providing conservation at the discrete 
level. 

We will now show that the ESWENO scheme Eq. (9) satisfies a broader definition of conservativeness 
which requires that not only the numerical flux itself, but also all moments of the flux against an arbitrary 
test function v(x,t) telescope across the domain. It can easily be shown that this property holds for the 
continuous problem. Indeed, multiplying Eq. (1) by an arbitrary test function v(x,t) that vanishes at the 
boundaries, and integrating with respect to x and t, one can obtain the weak form of the original differential 
equation 

l d t l 


f vuix 


J J ( uv t + fv x ) dxdr = 0. 


(57) 


0 lo 0 0 

To show that the ESWENO scheme possesses the same property at the discrete level, we multiply Eq. (9) 
by the discrete vector v T P (the discrete analog of integration in space), which yields 


v T Pu t + v t Qu + \ T SD\Vl = 0. 


(58) 


Integrating Eq. (58) with respect to time and taking into account that <5 and S are fully skew-symmetric 
and symmetric matrices, respectively, we have 


t t 

v T Pu|* — J (u T Pv( + u t Qv) dr + J u T D± SDiwdr = 0. (59) 

o o 

In the limit of infinite spatial resolution, the last term on the left hand side of Eq. (59) vanishes, and the 
ESWENO scheme approaches the continuous operator. Thus, the discrete P, Q, and S operators mimic the 
conservation property of the continuous operator, which provides that the discrete equation (59) satisfies a 
weak form similar to Eq. (57). 


V. ESWENO scheme for Systems of Equations 

The ESWENO scheme constructed for the one-dimensional constant coefficient wave equation (1) can be 
extended to hyperbolic systems. There are two major approaches to generalization of scalar finite difference 
WENO-type schemes to systems of equations. The first approach is to discretize each equation in the system 
by using the scalar WENO-type scheme. In this case, the WENO reconstruction is done for each component 
of the vector separately. The second technique is based on the characteristic decomposition. The main idea 
of this approach is to transform the undivided differences to the local characteristic fields and perform the 
scalar WENO reconstruction for each component of the vector of the characteristic variables. The numerical 
flux calculated this way is then transformed back into the physical space. We will show in this section that 
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the energy estimate can be obtained for the characteristic form of a hyperbolic system of constant coefficient 
equations, while it is not clear how to obtain a similar energy estimate when the system of equations is 
discretized in the component-wise fashion. 

Let us consider a one-dimensional hyperbolic system of n equations with constant coefficients 


u t + Au x = 0, (60) 

where ft is a vector of length n, A is a constant (n x n) matrix which has n real eigenvalues: Ai < . . . < A ra , 
and a complete set of linearly independent eigenvectors: rq, . . . , f n . By multiplying Eq. (60) by a matrix R 
composed of the eigenvectors of A , i.e. , R = (rq, . . . ,r n ), and by introducing a vector of the characteristic 
variables c = the system of equations (60) can be rewritten in the following diagonal form: 


c t + A c x = 0, (61) 

where A = diag(Ai, . . . , A„). The system of equations (61) is fully decoupled, so that each equation is a 
scalar constant coefficient wave equation of the form of Eq. (1). As in the scalar case, we assume that the 
systems of equations (60) and (61) are subject to periodic boundary conditions. 

In contrast to the scalar case where the wave velocity a is either positive or negative, the eigenvalues in 
the system of equations (61) may have different signs simultaneously. If this is the case, a flux splitting is 
used 

/ = /+ + r , 

where the matrix has only nonnegative eigenvalues, while the matrix has only negative eigen- 

values, where / = Au. For the characteristic form of the equations (61), the above splitting is simply 
equivalent to separating the equations with positive eigenvalues from those with negative eigenvalues. In the 
present analysis, the Lax-Friedrichs flux splitting given by Eq. (16) is used. 

Both Eqs. (60) and (61) can be approximated by using the ESWENO scheme. However, as will be shown 
next, an energy estimate can be obtained only for the characteristic form of the governing equations (61), 
while similar estimate is not available for Eq. (60). Discretizing each equation in Eq. (61) by using the 
scalar ESWENO scheme (9, 28, 32-34), the semi-discrete approximation of the characteristic form of the 
hyperbolic system of equations can be written as 


Ct. + 2A c C — — 2? a C, 


(62) 


where C = [c^, . . . , Cj \ . . . , c{ \ . . . , Cj"'] T , V c and P a are block matrices defined as follows: 


V, = 


Ai P-'Qi 

0 


0 

A nP^Qr. 




\iP~ 1 Df SiDi 0 


0 A n P~ 1 DjS n D 1 


(63) 


where P, Qi, and S), i = l,n are (J + 1) x (J + 1) matrices which are given by Eqs. (28) and (32-34), 
respectively. 

Hereafter in this section, without loss of generality, it is assumed that all eigenvalues of the matrix A are 
nonnegative, i.e., 

An > • • • > Ai > 0, (64) 

which is equivalent to using the flux splitting and considering only the positive flux f + = l A ^ A c, where 
|A| = diag(|Ai|, . . . , |A„|). For the flux /“, the ESWENO reconstruction and the corresponding matrices 
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Qi and S) can be obtained as mirror images of those used for f + with the respect to Xj + 1 . All the results 
presented in this section for the positive ESWENO flux f + remains valid for f~ as well. 

The smoothness indicators used for calculation of Qi and Si are functions of C; = [c^ Cj' l ] T . In other 
words, for tth component of the vector of the characteristic variables, the smoothness indicators, which in 
the scalar case are defined by Eq. (19), are given by 



whereas the same formulas (18) are used for calculation of the weight functions. Hence, 


(65) 


Qiy^Qi, Si ^ Si, if i^l, 

because the nonlinear operators Qi and Si depend on c , = [c^, . . . ,Cj^] T , and therefore are different for 
each component of the vector of the characteristic variables. 

We will now prove that the ESWENO scheme given by Eqs. (62-64) is stable in the energy sense. 

Theorem 2. The approximation (62-64) with the matrices Qi and Si, i = 1 ,n defined by Eqs. (28) and 
(32), respectively, is stable. 

Proof. An energy estimate for Eq. (61) can be obtained by employing the standard procedure used in 
the scalar case. Multiplying Eq. (62) with C T V leads to 

I±\\c\\ 2 v + C T QC = ~C T SC, ( 66 ) 

where V = diag(P, . . . , P ), ||C||)b = X)"=i ll c *l!p and 2 and S are the following block matrices: 



AiQi 

0 


' A rDfSrDr 

0 

Q = 

0 

^ nQn 

, 5 = 

0 

XnDfSnD! 


Adding Eq. (66) to its transpose, we have 

i\\Cf v + C t {Q + Q t )C = -C T {S + S T )C. 
at 

Taking into account that the matrices Q*, i = l,n are skew-symmetric and 


c? 

1 

0 


i 

i 

9 

0 

0 

nQn 


0 

^nQn 


( 68 ) 


(69) 


the second term on the left hand side of Eq. (68) vanishes, and Eq. (68) can be recast as follows: 


d_ 

dt 



\i{D 1 c i ) T (Si + Sj)D 1 C, < 0. 

i = 1 


(70) 


From Eq. (70), it follows that if Eq. (64) is valid, and each matrix Si, i = l,n is positive semidefinite, then 
the ESWENO scheme (62-64) is stable for systems of equations (61). □ 
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Remark 9. Note that the standard approach used for proving Theorem 2 does not provide a similar energy 
estimate when the system of equations (60) is discretized by the ESWENO scheme in a component by 
component fashion. Indeed, discretizing Eq. (60) by using the ESWENO scheme yields 

U t + V c u=-Au, (71) 


where U = , . . . , u'P , . . . , u ^ , . . . , u j n) ] r , V c and V a are defined as 


V c = 


V a = 


Oil P 1 Q 1 0,1 n P 1 Ql 

a n \P Q n . . . a nn P Q n 
a\\P~ l S\D\ ... a\ n P~ L D{ S\D\ 

a nl P~ 1 DjS n D 1 ... a nn P~ 1 DfS n D 1 


(72) 


Multiplying Eq. (71) by U T V and adding it to its transpose, we have 

i\Mv + U T (Q + Q T )U = - {ViUf (S + 5 t )2?iW, 

at 

where = diag(Z?i, . . . , D i), Q, Q T , S T , and S T are given by 



auQi 

ainQi 


auQi 

£ 

I 

Q = 

O'nlQn 

• • O'nnQn 

, Q T = 

O'lnQi 

• • Q"nnQ n 



anSi 

£ 

3 

Co 


auSf . 


s = 

^nl^n 

• • ®nn 

Co> 

II 

QinSf . 

• • ^nn 


(73) 


(74) 


From the above equation, it immediately follows that although Qi = —Qj, Vi = 1 ,n, the matrix Q is not 
skew-synnnetric, thus making the second term on the left hand side of Eq. (73) non- zero. Furthermore, the 
matrix S + S T is not positive semidefinite because of the presence of a, j . As a result, an energy estimate sim- 
ilar to one derived for the characteristic form of the hyperbolic system of equations (61) cannot be obtained 
when Eqs. (60) are approximated in the component-wise fashion. A conclusion that can be drawn from this 
analysis is that to guarantee the stability of the ESWENO scheme for hyperbolic systems, the character- 
istic decomposition should be used in the ESWENO reconstruction. As has been shown numerically in, 26 
the component- wise high-order WENO reconstruction gives rise to spurious oscillations for problems with 
strong discontinuities, while the characteristic- wise WENO reconstruction allows one to obtain essentially 
nonoscillatory solutions, which corroborates our theoretical results. 


VI. Numerical Results 

In this section, we will compare numerical results obtained with the third-order ESWENO and conven- 
tional WENO schemes. For all test problems considered, the CFL number is set to be 0.9, and the parameters 
k, m, e in Eqs. (18, 19) and 6 in Eq. (34) have been chosen based on the truncation error analysis presented 
in Section 4.3 

where J is the number of grid cells. The time derivative is approximated using the 3rd-order TVD 
Runge-Kutta method developed in. 25 
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k 

m 

e 6 

WENO 

2 

2 

1.0e-6 

ESWENO 

2 

1 

5.0e-3 J~ 2 


Table 1. The k , m, e, and 8 parameters for WENO and ESWENO schemes. 


VI. A. Scalar Linear Wave Equation 

We begin by verifying that the ESWENO scheme is third-order accurate for smooth problems. To check this 
property, we consider equation (1) with a = 0.5 and the following initial condition: 


«o(z) 


[0.5 + 0.5 cos(w(a: — x c ))] 4 , if \x — x c \ < a 
0, if \x — x c \ > <j 


(75) 


J 

Linear 

l|u-u ft || Loo 

Rate 

ESWENO 
l|u- u h \\ Loo 

Rate 

WENO 
II u- WlUoc 

Rate 

50 

3.98e-l 

- 

5.10e-l 

- 

6.82e-l 

- 

100 

1.73e-l 

1.20 

2.62e-l 

0.96 

5.02e-l 

0.44 

200 

3.99e-2 

2.11 

6.40e-2 

2.03 

2.72e-l 

0.88 

400 

5.78e-3 

2.78 

7.84e-3 

3.03 

9.43e-2 

1.53 

800 

7.41e-4 

2.97 

8.44e-4 

3.22 

2.88e-2 

1.71 

1600 

9.29e-5 

3.00 

9.67e-5 

3.13 

6.21e-3 

2.22 


Table 2. Loo error norms and their convergence rates obtained with the 3rd-order ESWENO, WENO, and linear 
schemes. 


J 

Linear 
1 u — W Li 

Rate 

ESWENO 
ll Q - u h |Ui 

Rate 

WENO 

l|u-u,,|| Ll 

Rate 

50 

6.60e-2 

- 

7.62e-2 

- 

1.08e-l 

- 

100 

2.49e-2 

1.41 

3.39e-2 

1.17 

6.74e-2 

0.68 

200 

5.05e-3 

2.30 

6.84e-3 

2.31 

2.71e-2 

1.31 

400 

7.03e-4 

2.85 

7.99e-4 

3.10 

6.02e-3 

2.17 

800 

8.93e-5 

2.98 

9.27e-5 

3.11 

1.18e-3 

2.35 

1600 

1.12e-5 

3.00 

1.13e-5 

3.03 

1.49e-4 

2.98 


Table 3. L i error norms and their convergence rates obtained with the 3rd-order ESWENO, WENO, and linear schemes. 


where u>, x c and a are equal to 57r, 0.5, and 0.2, respectively. The computational domain for this problem 
is set to be 0 < x < 1. Numerical solutions are calculated on a sequence of globally refined uniform grids 
and advanced in time up to t = 10 which corresponds to 5 periods in time. To verify the order property, the 
error convergence rate is evaluated as 


Rate = 


log ( || u - u Axi ||) — log ( || u - u A x2 ||) 
log(Aa:i) - log(Ax 2 ) 


(76) 


where u is the exact solution, u Aa:i and u AiC2 are the corresponding numerical solutions computed on grids 
with the grid spacings Aaq and Aa: 2 , respectively. 
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J 

llu-Uhllioo 

Rate 

l|u-u ft || Ll 

Rate 

50 

4.79e-l 

- 

7.23e-2 

- 

100 

2.31e-l 

1.05 

2.98e-2 

1.28 

200 

5.50e-2 

2.07 

5.82e-3 

2.35 

400 

7.38e-3 

2.90 

7.37e-4 

2.98 

800 

8.16e-4 

3.18 

9.01e-5 

3.03 

1600 

9.54e-5 

3.10 

1.12e-5 

3.01 


Table 4. Lqq and L\ error norms and their convergence rates obtained with the 3rd-order WENO scheme with k = 2, 

m = 1, e = 5 X 10 — 3 . 



Figure 1. Comparison of the 3rd-order linear, ESWENO, and conventional WENO schemes for the propagation of a 
cosine pulse on a uniform grid with J = 200 over 5 time periods. 


In regions where the solution is sufficiently smooth, both the WENO and ESWENO schemes convert to 
a 3rd-order linear, upwind-biased scheme with the following flux: 

fj + 1 = “g/i - 1 + g fj + 3 / 7 + 1 - (77) 

In the vicinity of strong discontinuities and unresolved features, WENO and ESWENO schemes introduce 
additional nonlinear artificial dissipation through the weight functions. The presence of the nonlinear dis- 
sipation terms make these schemes more dissipative than the target 3rd-order linear scheme (77), which 
provides the lower error bound that can be obtained with the WENO and ESWENO schemes for smooth 
solutions. The L ^ and L\ error norms and their convergence rates of the 3rd-order linear, ESWENO and 
WENO schemes are presented in Tables 2 and 3, respectively. In Table 4, we also present global grid re- 
finement results obtained with the WENO scheme whose tuning parameters have been set to be equal to 
those of the ESWENO scheme (k = 2, m = 1, e = 5 x 10 -3 ). As follows from these numerical results, the 
3rd-order ESWENO scheme significantly outperforms the conventional 3rd-order WENO scheme in terms of 
accuracy and is slightly more dissipative than the optimized WENO scheme because of the presence of the 
additional artificial dissipation term (43). The L i errors obtained with the 3rd-order linear, ESWENO, and 
optimized WENO schemes reach the design convergence rate, starting at J = 400 grid points, whereas the 
conventional WENO scheme demonstrates the design-order convergence rate only on the finest mesh with 
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J= 1600. 

For all grids considered, the maximum error occurs at the peak of the cosine function. As follows from 
Table 2, for the ESWENO scheme, the L ^ error convergence rate is close to three on sufficiently fine meshes, 
thus indicating that the ESWENO scheme is 3rd-order accurate at the smooth extrema. However, the L 0 0 
norm of the error obtained with the conventional WENO scheme exhibits only the 2nd-order convergence 
rate at the finest mesh and is almost two orders of magnitude larger than that of the ESWENO scheme at 
J = 1600. The accuracy of the WENO scheme can be drastically improved and the design order can be 
recovered on moderate grids if the tuning parameters are set to be equal to their optimized values, as one can 
see in Table 4. The comparison of solutions obtained with the linear, ESWENO, and conventional WENO 
scheme on a 201-point grid are shown in Fig. 1. The solution obtained with the optimized WENO scheme 
is practically indistinguishable from that of the ESWENO scheme and therefore is not presented in Fig. 1. 
Figure 1 shows that the cosine pulse amplitude computed using the conventional 3rd-order WENO scheme 
is much lower than that of the exact solution and the amplitude obtained with the 3rd-order ESWENO 
scheme, thus showing that the ESWENO scheme is much less dissipative than its conventional counterpart. 
Another observation is that for smooth solutions, the ESWENO scheme provides almost the same accuracy 
as the target, linear 3rd-order scheme. 



Time 


Figure 2. Time histories of the rightmost eigenvalue of the symmetric part of the ESWENO and WENO operators. 

The next question that we would like to address concerns the stability of WENO and ESWENO schemes. 
To the authors’ knowledge, no proof has been reported on the stability of WENO schemes for problems with 
discontinuous solutions. As has been discussed in Section 4, for the linear convection equation with piecewise 
continuous initial conditions, the ESWENO scheme is stable in the energy sense. A sufficient condition for 
the stability is that —^(Des + £>es) is negative semidefinite, where Des is the ESWENO discrete operator. 
This property implies that all eigenvalues of the symmetric part of the ESWENO operator are nonpositive. 
It should be noted that the symmetric part of the WENO operator may have positive eigenvalues as has 
been discussed in Section 3.2. To verify these properties, time histories of the rightmost eigenvalue of the 
symmetric part of the ESWENO and WENO operators computed on a 101-point grid are compared in Fig. 
2. On the entire time interval considered, the maximum eigenvalue of the -1 (Des + ^es) matrix is equal 
to zero to the order of the round-off error. In contrast to the ESWENO scheme, the symmetric part of the 
conventional 3rd-order WENO operator has positive eigenvalues of the order of 0(1), as one can see in Fig. 
2. This is due to the fact that the symmetric part of the WENO operator is not negative semidefinite, if 
there are unresolved features in the computational domain. The unstable modes of the WENO operator 
can be significantly suppressed by adjusting to and e in Eqs. (18, 19). When the tuning parameters are 
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equal to their optimized values {in = 1, e = 5 x 10” 3 ), the maximum positive eigenvalue of the symmetric 
part of the WENO operator is shifted toward the imaginary axis and monotonically decreases from 10 ~ 2 to 
10~ 5 , which indicates that this choice of the turning parameters does not only improve the accuracy but also 
stabilizes the scheme. Another important observation is that the rightmost eigenvalue of the conventional 

ESWENO 



Figure 3. Time histories of the WENO and ESWENO weight functions, a; 0 , at the cosine peak. 


WENO operator exhibits highly oscillatory behavior. These liigh-amplitude oscillations indicate that the 
presence of eigenvalues in the right-half plane makes the conventional WENO scheme locally unstable. As 
a result, the WENO scheme biases the stencil away form the local instabilities caused by these eigenvalues, 
thus introducing additional dissipation which shifts the positive eigenvalues closer to the imaginary axis. 
One can see this behavior of the WENO weight functions in Fig. (3) which shows time histories of the 
WENO and ESWENO weights, w°, at the cosine peak. As the cosine peak, which is not fully resolved on 
the 101-point mesh, advances forward in the computational domain, the eigenvalues corresponding to those 
grid points that are in the vicinity of the unresolved feature are shifted to the right. This pattern persists 
over the entire time interval considered. In contrast to the WENO scheme, the ESWENO weight function, 
w°, slightly oscillates and approaches to 2/3 which is its optimal value for smooth solutions. Note that the 
presence of the positive eigenvalues does not imply that the 3rd-order WENO scheme is globally unstable, 
because these eigenvalues correspond to different grid points at different moments of time. It should also be 
noted that the results obtained with the optimized WENO scheme are very similar to those of the ESWENO 
scheme and therefore are not presented hereafter. 

The second test problem has been chosen to check whether the ESWENO scheme provides essentially 
nonoscillatory solutions for problems with strong discontinuities. The same linear convection equation (1) 
with a = 0.5 is used as a test problem. However, in contrast to the previous test case, the initial condition 
is set to be a piecewise continuous function given by 


uq{x) 


1, if \x — x c \ < <j 
0, otherwise 


(78) 


where x c = 0.5 and a = 0.2. Numerical solutions obtained with the 3rd-order linear, ESWENO, and WENO 
schemes on a uniform grid with 101 points at time t = 2 are compared with the exact solution in Fig. 4. As 
expected, the solution computed using the linear upwind-biased scheme is oscillatory and has large-amplitude 
overshoots and undershoots near the discontinuities. The WENO results are free of spurious oscillations and 
the most dissipative among those presented in Fig. 4. The ESWENO solution is less dissipative than that 
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obtained with the conventional WENO scheme and essentially nonoscillatory near the strong discontinuities. 



x 


Figure 4. Comparison of the 3rd-order linear, ESWENO, and WENO schemes for the propagation of a square pulse 
on a uniform grid with J = 100 over 1 time period. 


Similar to the previous test case, we numerically analyze spectra of the symmetric part of the WENO 
and ESWENO finite difference operators for the linear convection equation (1) with the discontinuous initial 
condition (58). As has been shown in Sections 3.2 and 4.1, the symmetric part of the WENO operator can 
have positive eigenvalues, whereas the spectrum of the symmetric part of the ESWENO operator is always 
located in the left-half plane regardless of whether the solution is continuous or discontinuous. To verify these 
properties, time histories of the rightmost eigenvalue of the symmetric part of the ESWENO and WENO 
operators are depicted in Fig. 5. As expected, the maximum eigenvalue of — ^(-Des + D^ s ) is of the order of 
the machine zero. The symmetric part of the conventional WENO operator has positive eigenvalues at each 
time step. The largest positive eigenvalue of — |(Dw + D oscillates at the beginning and then decreases 
reaching its minimum value at t « 5, thus indicating that the WENO dissipation operator suppresses the 
unstable modes. For later times, the largest positive eigenvalue gradually increases up to 0.05 and stays at 
this level until the end of the time interval considered. It should be noted that for this test problem the 
conventional WENO scheme is stable in the sense that the numerical solution remains bounded and does 
not grow exponentially. 

VI. B. The 1-D Euler Equations 

As has been proven in Section 5, the 3rd-order ESWENO scheme is energy stable for a system of linear 
hyperbolic equations with periodic boundary conditions. Although, there is no proof that the ESWENO 
scheme is stable for nonlinear conservation laws with nonperiodic boundary conditions, it is interesting to 
see how the scheme performs in this case. In connection with this, we consider the quasi- 1-D Euler equations 
for a perfect gas, which are given by 


3U i 
dt "T" 

H | ►Tj 

II 

J 

G 





P 


pu 


pu 

u = 

pu 

E 

, F = 

pu 2 +p 
(E +p)u 

0 

II 

pu 2 

(E +p)u 


where A = A(x) is the cross-sectional area of a quasi-l-D nozzle. We will first use the ESWENO and WENO 
schemes to solve steady subsonic and transonic quasi-l-D nozzle flows and then compare solutions obtained 
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Figure 5. 


Time histories of the rightmost eigenvalue of the symmetric part of the ESWENO and WENO operators. 


with these schemes for an essentially time-dependent problem. For the steady problems, both the ESWENO 
and WENO residuals were driven below 10 -11 . 

In Section 5, it has been shown that to preserve the stability of the ESWENO scheme for systems of 
hyperbolic conservation laws, the ESWENO reconstruction should be carried out in local characteristic 
fields, which has been implemented as recommended in. 27 First, the Roe average at j + \ is used to form the 


. i composed of the left and right eigenvectors of 9F/9U. Then, all quantities required 


matrices Rj+i and R 

to compute the ESWENO flux are multiplied by J?' 1 , . 

+ 2 

flux spitting: 


After that we perform the following Lax-Friedrichs 


— 2 (F A A max U) , A max — diag (A max , A^ ax , Aj^ ax ) 
Am ax = max A" i , n = U3 

0 <J<J •*' 2 


(80) 


and use the scalar ESWENO reconstruction for each component of the characteristics variables. Once the 
positive and negative ESWENO fluxes at j + \ have been formed, they are transformed back into the physical 
space by left multiplying them with Rj + i ■ Note that the same local characteristic decomposition is used for 
the WENO scheme, but instead of the ESWENO reconstruction, the WENO reconstruction is used. 

The steady state isentropic flow through a quasi-l-D nozzle with the cross-sectional area 


A(x) = 1 — 0.8a:(l — x), 0 < x < 1 


has been chosen as a test problem. The inflow Mach number is set to be 0.5 and the pressure at x = 1 is 
assumed to be equal to that at x = 0. Under these conditions, the flow inside the nozzle is fully subsonic. 
Global grid refinement studies for the linear, WENO, and ESWENO schemes are presented in Tables 5 
and 6. As follows from this comparison, the L ^ and L\ error norms obtained with the ESWENO scheme 
are about an order of magnitude less than those of the WENO scheme and very close to their theoretical 
limits obtained with the 3rd-order linear upwind-biased scheme. Note that the WENO scheme exhibits the 
4th-order convergence rate on sufficiently fine meshes, which is due to the fact that on coarse meshes, the 
conditions (37, 40) are not met, thus making the WENO scheme too dissipative. On finer meshes these 
conditions are satisfied, and the error norm rapidly approaches its theoretical limit at the rate which is 
greater than its asymptotic limit. 
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J 

Linear 
||u - wlkoo 

Rate 

ESWENO 
l|u— u h \\ Loo 

Rate 

WENO 

llu-Uhlk*. 

Rate 

50 

6.86e-6 

- 

1.57e-5 

- 

6.34e-4 

- 

100 

8.43e-7 

3.02 

1.55e-6 

3.35 

9.67e-5 

2.71 

200 

1.04e-7 

2.21 

1.79e-7 

3.12 

6.29e-6 

3.94 

400 

1.296-8 

2.82 

2.18e-8 

3.04 

3.81e-7 

4.05 

800 

1.61e-9 

2.97 

2.70e-9 

3.01 

2.08e-8 

4.19 


Table 5. Loo error norms and their convergence rates obtained with the 3rd-order linear, ESWENO and WENO schemes 
for the subsonic quasi-l-D nozzle problem. 


N 

Linear 

llu-Uhlk, 

Rate 

ESWENO 

lla-ufclk 

Rate 

WENO 

llu-u^k 

Rate 

50 

2.73e-6 

- 

5.97e-6 

- 

9.37e-5 

- 

100 

3.38e-7 

3.02 

5.97e-7 

3.33 

1.25e-5 

2.90 

200 

4.19e-8 

3.01 

6.93e-8 

3.10 

1.23e-6 

3.34 

400 

5.21e-9 

3.01 

8.48e-9 

3.03 

1. Ole-7 

3.60 

800 

6.50e-10 

3.00 

1.05e-9 

3.01 

6.86e-9 

3.89 


Table 6. L i error norms and their convergence rates obtained with the 3rd-order linear, ESWENO and WENO schemes 
for the subsonic quasi-l-D nozzle problem. 


The next test problem is the steady transonic flow through a quasi-l-D nozzle with the following cross- 
sectional area: 

A{x) = 1.398 + 0.347 tanh(0.8x — 4), 0 < x < 10. 

The Mach number at x = 0 is 1.5, and the outflow conditions have been chosen so that the shock is 
located at x = 5. Density profiles calculated using the 3rd-order ESWENO, WENO, and linear schemes 
on a 51-point grid and the exact solution are shown in Fig. 6. Similar to the square pulse test problem, 
the numerical solution obtained with the 3rd-order linear scheme oscillates near the shock. For both the 
ESWENO and WENO schemes, the captured shock is smeared over three grid cells, and the numerical 
solutions are essentially nonoscillatory and agree very well with the exact solution. To quantify the accuracy 
of each scheme in regions where the solution is smooth, we measure the L\ error norm downstream of the 
shock at x > 6, which is presented in Table 7. As in the previous test problems considered, the ESWENO 
scheme is about an order of magnitude more accurate than the conventional WENO scheme on coarse and 
moderate grids. Note that on the coarsest 51-point mesh, the ESWENO scheme is more accurate than the 
target linear scheme, which is due to the presence of Gibbs oscillations in the solution obtained with the 
linear scheme. It is not surprising that all three schemes provide the design order of accuracy, because as has 
been shown in, 28 for the steady quasi-ID nozzle problem the first-order error component generated by the 
shock-capturing procedure is localized near the shock, which is not the case for time-dependent and steady 
multidimensional problems. 29 

The last test problem considered is the interaction of a moving shock with smooth density fluctuations. 
The solution of this benchmark problem contains both strong discontinuities and smooth structures, which 
is very well suited for testing high-order shock-capturing schemes. The governing equations are the time- 
dependent 1-D Euler equations (79) with G = 0, which are subject to the following initial conditions: 

(p,u,p) =(3.857134,2.629369,10.33333), for x < -4 , . 

(p,u,p) = (1 + 0.2 sin 5a;, 0, 1) , for x > — 4. 

The governing equations are integrated in time up to t = 1.8. The exact solution to this problem is not 
available. Therefore, a numerical solution obtained with the conventional 3rd-order WENO scheme on a 
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Figure 6. Comparison of the 3rd-order ESWENO, WENO, and linear schemes for the steady transonic flow through 
the quasi-l-D nozzle. 


J 

Linear 

l|u~ WlUx 

Rate 

ESWENO 
l|u- wlk 

Rate 

WENO 

llu-Uhlli,! 

Rate 

50 

2. Ole-5 

- 

1.51e-5 

- 

1.33e-4 

- 

100 

9.69e-7 

4.38 

1.40e-6 

3.43 

1.56e-5 

2.90 

200 

1.17e-7 

3.04 

1.67e-7 

3.06 

1.53e-6 

3.34 

400 

1.44e-8 

3.03 

2.05e-8 

3.03 

1.13e-7 

3.60 

800 

1.78e-9 

3.01 

2.53e-9 

3.02 

6.86e-9 

3.89 


Table 7. L i error norms and their convergence rates measured downstream of the shock at x > 6, which were computed 
using the linear, ESWENO and WENO schemes for the transonic quasi-l-D nozzle problem. 


uniform grid with J = 10000 grid cells is used as a reference solution. 

The parameters k, m, e, and <5 for the ESWENO and WENO schemes have been set to the same values 
used in all previous test problems, which are given in Table 1. The local characteristic decomposition based 
on the global Lax-Friedrichs flux splitting (80) is used for both the ESWENO and WENO reconstructions. 
Numerical solutions obtained with the 3rd-order ESWENO and WENO schemes on a 801-point grid at 
t = 1.8 are compared with the “exact” reference solution in Fig. 7. Note that both the ESWENO and 
WENO solutions are free of spurious oscillations. Figure 7 shows that the 3rd-order ESWENO scheme 
is more accurate than its conventional counterpart and provides much better resolution in a region right 
upstream of the moving shock. As has been discussed above, two major factors that make the WENO 
scheme too dissipative are the choice of the parameters k, m, and e in Eqs. (18, 19) and the presence of 
positive eigenvalues in the spectrum of the symmetric part of the WENO operator. Note that for this test 
problem, the 3rd-order target, linear scheme is unstable at CFL=0.9. 

VII. Conclusions 

A new third-order ESWENO finite difference scheme is developed for scalar and vector linear hyperbolic 
equations with discontinuous solutions. The new scheme combines the standard third-order WENO scheme 
developed by Jiang and Shu and an additional nonlinear artificial dissipation term. The additional term is 
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Rho 




Figure 7. Density profiles computed with the 3-rd-order ESWENO (left) and WENO schemes on a uniform grid with 
801 points for the density sine wave/shock interaction problem. 


third-order accurate for smooth solutions including smooth extrema. This term guarantees that the new 
scheme is stable in the /^-energy norm for both continuous and discontinuous solutions, whereas for the 
conventional WENO scheme, an energy estimate is not available. The distinctive feature of the ESWENO 
scheme is that the stability has been proven for any form of weight functions that satisfy Eq. (20). For 
the same choice of the weight function tuning parameters, the ESWENO scheme is slightly more dissipative 
than the WENO scheme, because of the presence of the additional dissipation term. Based on the rigorous 
truncation error analysis, the weight functions of the new scheme have been optimized, which drastically 
improves the accuracy of the ESWENO scheme as compared to its conventional counterpart, while essentially 
preserving the nonoscillatory properties near strong discontinuities. The spectrum of the symmetric part of 
the ESWENO operator is always located in the left half-plane, whereas the symmetric part of the WENO 
operator has positive eigenvalues, thus indicating that the WENO scheme may become locally unstable. 
The numerical experiments have shown that the ESWENO scheme with the optimized weights is much 
more accurate than the conventional WENO scheme in regions where the solution is smooth and provides 
essentially nonoscillatory solutions near strong discontinuities. 
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